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Abstract 

We present the first fully non-linear calculation of inflaton decay. We map 
inflaton decay onto an equivalent classical problem and solve the latter nu- 
merically. In the Xcj) 4 model, we find that parametric resonance develops 
slower and ends at smaller values of fluctuating fields, as compared to esti- 
mates existing in literature. We also observe a number of qualitatively new 
phenomena, including a stage of semiclassical thermalization, during which 
the decay of inflaton is essentially as effective as during the resonance stage. 
PACS numbers: 98.80. Cq, 14.80.Mz, 05.70.Fh 
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Amplification of quantum fluctuations and transition to semiclassical behavior are famil- 
iar subjects in inflationary cosmology (for a review of the latter, see Ref. 0). Semiclassical 
fluctuations produced during inflation can explain the observed large scale structure of the 
Universe [[J. The theory of their production is by now well developed; for a recent rather 
general presentation, see Ref. 0. After inflation ends, the scalar field (inflaton) that was 
driving it begins to oscillate. These oscillations are thought to lead to particle production, 
inflaton decay, and eventually to reheating of the Universe. 

It was only recently realized that phenomena related to large occupation numbers of Bose 
particles produced by inflaton decay, such as parametric resonance, can be significant in the 
reheating process That means that semiclassical phenomena may play as important 

role in reheating as they do in the structure formation. Of course, the eventual thermalized 
state will include modes with small occupation numbers, which are essentially quantal, but 
initial stages of the reheating process should admit a semiclassical description. 

However, semiclassical description of fluctuations produced by inflaton decay has not 
been obtained. The existing treatments of parametric amplification employ the stan- 
dard methods of analyzing particle production, based on Bogoliubov transformation. These 
methods allow one to take into account some of the back reaction effects of produced par- 
ticles on the oscillating zero-momentum (k = 0) mode, such as time-dependence of the 
frequency of oscillations, but cannot reliably take into account many other non-linear effects 
that become important certainly no later than the frequency change. These effects include 
scattering of produced particles off the k = mode, which knocks particles out of the k = 
state; rescattering of the produced particles, which populates non-resonant modes as well as 
returns particles back to k = (Bose condensation) 0, etc. We generically refer to these 
effects as rescattering; we will see that they begin as essentially classical interactions. 

Order of magnitude estimates of the amplitude squared of produced fluctuations (<50) 2 
at the end of the resonance stage were obtained in Refs. {§-[(|. It is important to know {5(f)) 2 
more precisely, because in realistic inflationary models, its magnitude determines whether 
any symmetries are restored in non-thermal regime after the resonance stage ||. This, in 
turn, requires taking into account all rescattering effects. After the parametric resonance 
ends, rescattering leads to chaotic behavior and begins to thermalize the system, while it 
is still in the semiclassical regime. This semiclassical thermalization, which should proceed 
rather fast, due to the still large occupation numbers, is very important for determining 
the reheating temperature. However, it is completely absent from the standard treatment 
as well as from the one augmented by time-dependent Hartree approximation |J (in 
the Hartree approximation, chaotic behavior characteristic of a non-linear classical system 
is lost). 
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It is clear, on the other hand, that if one succeeds in finding an adequate semiclassical 
description of fluctuations produced by inflaton decay, the decay can be studied via the 
classical equations of motion and thus treated as a fully non-linear problem, without any 
further assumptions. That classical problem can be turned over to a computer. To obtain 
the semiclassical description, one has to find the initial conditions for the classical equations 
of motion, corresponding to the original quantum state. We find these initial conditions 
using the following observation. Because inflationary models have small coupling constants, 
fluctuating modes can grow to large occupation numbers before their interactions become 
important. In other words, the transition to classical behavior occurs in a linear regime with 
respect to fluctuations. We can then use the linear theory of quantum-to-classical transitions 
H to obtain the initial conditions for subsequent non-linear classical evolution. Possibility 
of a semiclassical formulation of the problem of inflaton decay was noted in a recent paper 
by Son || ; however, the correct initial conditions for the classical problem were not obtained 
there. 

In this Letter we report the first results of our numerical study of the semiclassical decay 
of inflaton, for the simplest model of a real scalar field with the A0 4 potential. Initial condi- 
tions were chosen to correspond to conformal vacuum for fluctuations at the end of inflation 
|T0| . Numbers for various quantities of physical interest are presented below; because ours 
is the first fully non-linear treatment of the problem, they differ from estimates existing in 
the literature. We also observe a number of qualitatively new phenomena, including the 
stage of inflaton decay when the resonance peaks in power spectrum begin to interact and 
smear out; we call this stage semiclassical thermalization. We expect this stage on general 
grounds, as noted above, but it was not observed previously. We find that during semiclassi- 
cal thermalization inflaton decay is essentially as effective as during the stage of parametric 
resonance. 

We consider a real scalar field minimally coupled to gravity, in a Friedmann-Robertson- 
Walker universe. We work with rescaled variables: time r, such that a(r)dr = 
V\(j)o(0)a(0)dt, coordinate £ = \/A0o(O)a(O)x, an d field ip = 0a(r)/0 o (O)a(O), where a(r) 
is the scale factor, and </>o(0) is the zero- momentum mode at t — r — 0. In this Letter we 
consider only the massless inflaton case. (The method itself applies to other models, as well 
as to various laboratory systems.) The action is 

Q 1 IaHA \ l (■ M 2 ( V ^) 2 <f 

where dots denote derivatives with respect to r. The parameter a/ a = h is the rescaled 
Hubble parameter, h{r) = H (t) / \^X(f) (0) . For the massless inflaton, H(0) ~ -\/A0q(O) / M P1 , 
so that h(0) is independent of A and is of order one. The equation of motion following from 



(1) 
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(I) is 

if - V\<p - % + = . (2) 

When the massless field begins to oscillate, the Universe becomes radiation dominated, so 
that a = 0, and the corresponding term drops out of Eq. (0). Although the coupling A 
does not appear in (0), it appears in the action and hence regulates the magnitude of initial 
quantum fluctuations in non-zero modes relative to the magnitude of the zero mode. We 
have checked numerically that in the massless case the condition for a semiclassical regime 
to occur is precisely the smallness of the coupling A. The zero mode is semiclassical from 
the beginning. Its initial value is <fo(0) = 1, by virtue of our rescaling. It is possible to 
choose <fo(0) = as a definition of the moment of time when the oscillations start; we do 
that. We make no further assumptions about the zero mode: its time dependence will come 
out of the solution to the full non-linear problem, once we specify the initial conditions for 
non-zero modes. 

According to the preceding discussion, to obtain the initial data for the classical problem, 
we can linearize Eq. (||) with respect to fluctuations. The linearized equation of motion for 
the Fourier transform of the fluctuation field (k ^ 0) is 

^k + ^fe( r Vk = ( 3 ) 

where u>l(r) = 3<Po( T ) + k 2 . The operator solution to Eq. @ is v 9 k( r ) = fk( T )b\ [ (0) + 
fk( T W wnere c-number function /^(t) satisfies Eq. (|3|) and the initial conditions 

A(0) = (2^0)) 2 ! fk(0) = HM0) + h(0)}f k (0). (4) 

Operators and b are the creation and annihilation operators defined at zero time and nor- 
malized according to [6^, 6^,] = 5(k— k'). A mode becomes semiclassical if the corresponding 
function fkij) grows (exponentially) for large r (and so does, then, the corresponding occu- 
pation number). In that case, it can be shown that up to exponentially small corrections, 
fk at large r can be made real by a time-independent phase rotation. When the exponen- 
tially small terms are neglected, v ? i c ( r ) an d its canonical momentum commute, and therefore 
V 9 k( r ) a ^ l &r S e times can be regarded as a random classical variable. This is the semiclassical 
description. 

The distribution of values of </?k( r ) is obtained from the solution to the linearized quan- 
tum problem. It depends on the quantum state for fluctuations existing at the end of 
inflation. Modes with k ~ 1 (in our rescaled units), which we are interested in, could only 
be produced towards the end of the de Sitter epoch, so the quantum state for these modes 
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should not be much different from conformal vacuum. Below we present results obtained 



using conformal vacuum as the original quantum state |T(J. The corresponding classical 
distribution function at large times is, cf. Ref. ||, 

ffy, p; r] = N exp (-j'^kLd 3 ^ 5 (f k (r)<p k - f k (r)p k ) (5) 

where the delta function is a shorthand for the product of functional delta functions for real 
and imaginary parts, and the prime on the integral reminds us that the integral is taken 
over half of the values of k (since these give all independent (p k ) and does not include k = 0; 
TV is a time-independent normalization. A useful point here is that (|5[) satisfies the classical 
Liouville equation also at small times. So, we can use it as an initial condition even at r = 0, 
where and its derivative are given by Eq. 

We thus use a discretized version of the distribution 

V[p(0)] = Nexp (-^J ^(0)|^ k (0)| 2 d 3 ^) (6) 

to generate random initial values of ip k . Once a value of p k is generated, the corre- 
sponding "velocity" is determined uniquely, according to eqs.(^) and (||), as <^ k (0) = 
[—iu)k(0) + h(0)] p k (0). Together with the initial conditions for the zero mode, this forms 
(after an inverse Fourier transform) an initial data problem for the fully non-linear equation 

(0)- 

We have integrated numerically Eq. (0) in a box of finite size L with the above initial 
conditions. We have varied the size of the box and the number of grid points iV to make 
sure that we are close to the continuum limit. We have varied the coupling constant A over 
the range 10 -10 -1. Here we present the results for L = 16n, N = 128 3 and A = 10~ 4 . All 
quantities we measured are averages of some sort; for these quantities, there is no need in 
further averaging over initial conditions. 

One important quantity that can be measured is the power spectrum of fluctuations, 
P(k). It is proportional to ( P k ^P\ s _ avera S e d over the direction of k and is normalized in such 
a way that Parseval's theorem reads 

£kP(k) = J- / d^[p(o - p f . (?) 



L 3 

Thus normalized power spectrum does not depend on the size of the (large) box. It is 
presented in Fig. |I] for several moments of time. Resonance peaks develop at r < 200 in 
the following order: first at k ~ 1, somewhat later at k close to zero, and still later at larger 
momenta. For comparison, we have run the problem linearized with respect to fluctuations. 
Only the initial development of the first peak was unaffected by linearization. The peak 
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at k ps 2 was barely visible, the others were not seen at all. We conclude that all peaks 
in Fig. |l| except the one £it k ~ 1 are due to rescattering. The peak at k close to zero 
is not expected from the approximate linear analysis based on the Mathieu equation. We 
interpret it as a result of the rescattering process in which a particle from the first peak 
transfers some of its momentum to a particle from the k = condensate |H |. The positions 
of the k ps 2 and higher peaks are close to those predicted by the Mathieu equation, but 
their widths and magnitudes are not. We interpret them as another result of rescattering: 
for example, two particles from the first peak, with k ps 1, together with two particles from 
the condensate produce two particles with k ps 2. (That process may be more efficient than 
the 6^2 process, which uses only condensate particles and is approximately described by 
the Mathieu equation.) Thus, rescattering processes start to play some role already at the 
parametric resonance stage. 

Another interesting quantity is the integral of the power spectrum of produced fluctua- 
tions over all (non-zero) k. This is simply the variance Var(<^) = ((<p(x) — (y?)) 2 ) and can be 
measured independently. The brackets denote averaging over the spatial lattice, so (0) is just 
the zero-momentum mode. The lattice effectively subtracts contribution of high-momentum 
modes, so at small A our Var(y?) approximates an already subtracted continuum quantity, 
equal to the continuum Var(<^) minus its value at zero time. The behavior of Var(t^) with 
time is shown in Fig. 0. Initially it grows exponentially; this is the stage of parametric 
resonance. Using logarithmic plots, we have found the following analytical fit for this stage 

In Var(v?) = 2 fir + In A - 8. (8) 

Because Var(yj) oscillates, the choice of constant here is somewhat arbitrary; we chose it in 
such a way that the line (|8|) passes through the maximum of the amplitude of Var(yj) at 
the end of the resonance stage (located at r = t* ps 200 if A = 10~ 4 ). For the rate of the 
exponential growth 2/i we obtain 2\i = 0.07, for all small enough A. This value is smaller 
than the one obtained from the Mathieu equation arising when the time dependence of the 
zero mode (an elliptic function in this regime) is replaced by sine with the same period. 
When we make this replacement in our program, we obtain for the rate the value of 0.11, in 
our time units, as expected. As another test, we have varied the size of the box and verified 
that we do not miss the fastest growing mode. 

The value of Var(^) in the maximum at the end of the resonance stage is approximately 
0.07 and does not depend on A when A 1. Thus, for general (small) A the resonance ends 
at time r* = 76. — 14.3 In A. Extrapolating to a realistic value A = 10~ 13 , we find ps 500. 
To obtain the variance of the physical field at time r, we need, at large r, to multiply 
Var(yj) by 3Mpj/27rr 2 . Hence, the maximum variance of 0, at the end of the resonance stage, 
is 
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Var(0) = I. x 1(T 7 M^ 



(9) 



For the effective temperature T e g, defined by equating (§) to T e 2 ff /12, this gives T e g = 1. x 
10 _3 Mpi. Notice that T e $ does not depend on the initial value of the inflaton field after 
inflation. The smallness of T e g is a result of two factors not taken into account previously: 
Var(<£>) <C 1 and very large r*. In realistic models including more fields, T e ff determines if 
any symmetries are restored in the non-thermal regime after the resonance stage ||. If in 
those models T e g- remains three orders of magnitude lower than the Planck scale, then, for 
example, non-thermal restoration of GUT symmetry may be prevented. On the other hand, 
we see that T e g is much larger than the reheating temperature that would be obtained if one 
neglected the effect of large occupation numbers, which confirms the main point of Refs. M . 

Finally, the time dependence of the zero-momentum mode is shown in Fig. [| We can 
see that after decreasing during the parametric resonance, it rebounces twice during time 
t = 200 to 300. We attribute this to Bose condensation caused by rescattering. Comparison 
of Figs. |2| and [3] suggests that it is Bose-condensation that terminates parametric resonance. 
However, the magnitude of the zero mode continues to decrease significantly, approximately 
as r -1 / 3 (as opposed to a slower t 1 / 12 decay found in ref. Returning to Fig. 

|l] we see that the peaks smear out at this stage. This is the onset of classical chaos or, in 
the language of particle physics, stimulated scattering, decay and annihilation due to large 
occupation numbers. This process is rather effective. In our model the energy scales as 
the forth power of the field. So, for A = 10~ 4 , at time r m 400 about 70% of energy is 
in fluctuations. For smaller A, the onset of this chaotic regime, which we call semiclassical 
thermalization, is delayed, but its rate with respect to time r does not decrease. The reason 
is that the strength of rescattering is determined not by A itself but by the product of A and 
a relevant occupation number. 

In conclusion, we have done the first fully non-linear calculation of inflaton decay. We 
have mapped inflaton decay onto an equivalent classical problem and solved the classical 
problem numerically. In the A0 4 model, we have found that parametric resonance develops 
slower and ends at smaller values of fluctuating fields, as compared to estimates existing in 
the literature. We have also observed a number of qualitatively new phenomena, including a 
stage of semiclassical thermalization (chaos), during which the decay of inflaton is essentially 
as effective as during the resonance stage. 
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Science Foundation under grant PHY 95-01458, and by the Alfred P. Sloan Foundation. The 
work of I.T. was supported by DOE grant DE-AC02-76ER01545 at Ohio State. 
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FIG. 1. Power spectrum of fluctuations at successive moments of time. 
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FIG. 2. Variance of the scalar field as a function of time. 
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FIG. 3. Time dependence of the zero-momentum mode. 
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